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Abstract 

Based on the standard many-fermion field theory, we construct models describing ultracold 
fermions in a ID optical lattice by implementing a mode expansion of the fermionic field operator 
where modes, in addition to space localisation, take into account the quantum numbers inher- 
ent in local fermion interactions. The resulting models are generalised Hubbard Hamiltonians 
whose interaction parameters are derived by a fully-analytical calculation. The special interest 
for this derivation resides in its model-generating capability and in the flexibility of the trapping 
techniques that allow the tuning of the Hamiltonian interaction parameters over a wide range of 
values. While the Hubbard Hamiltonian is recovered in a very low-density regime, in general, 
far more complicated Hamiltonians characterise high-density regimes, revealing a rich scenario for 
both the phenomenology of interacting trapped fermions and the experimental realization of de- 
vices for quantum information processing. As a first example of the different situations that may 
arise beyond the models well known in the literature (the unpolarised-spin fermion model and the 
noninteracting spin-polarised fermion model), we derive a Rotational Hubbard Hamiltonian de- 
scribing the local rotational activity of spin-polarised fermions. Based on a standard techniques we 
obtain the mean-field version of our model Hamiltonian and show how different dynamical algebras 
characterize the case of attractive and repulsive two-body potentials. 

PACS numbers: 71.10.Fd, 05.30.Fk, 03.75.Ss 
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I. INTRODUCTION 

Since the Bose-Einstein condensation of alkali atoms in magnetic traps , a massive 
experimental and theoretical effort has been dedicated to the investigation of confined atoms 
in the extremely low-temperature regime (for a review see 

The flexibility of optical trapping techniques has suggest ed the devise of different config- 
urations (lattices Q, , 0, [l^ 5 superlattices, etc. 12, 13, 14, 15, opening a vast 
scenario of research. The ability to tune atomic interactions via a magnetic field (Feshbach 
resonance [17[), along with the proposal of single atom trap loading techniques [18[, has 
proven to be of capital importance for ultracold fermions physics, yielding the possibility to 
study fundamental aspects of superfluidity (BCS-BEC crossover, see e. 
and envisaging new perspectives in quantum information processing [2^ IsJ 

The present work focuses on the theoretical investigation of the properties of (few) 
fermionic ultracold atoms loaded into a ID optical lattice, where global confinement is 
ensured by a magnetic trap. The description of such a physical system can be naturally 
performed in terms of a generalised Hubbard Hamiltonian (gHH) which is deduced from a 
general field-theoretic Hamiltonian with two body interaction 25[. At this stage, particular 
care must be taken in the choice of the function basis for the field operator expansion. Al- 
though the symmetries of the system can provide selection rules that reduce the involvement 
of the gHH, the resulting coefficient structure is very rich and, as a direct consequence, the 
Hamiltonian hardly tractable. Nevertheless, the generality of the model gives rise to a wealth 
of sub-models, depending upon different approximations and regimes. The guideline to find 
simplified Hamiltonians is given by the thorough analysis of the gHH coefficient structure. 

From this perspective the analytical knowledge of the coefficients is a powerful tool to 
establish the physical relevance of different sub-models in the various situations that may be 
conceived in the framework of the trapped ultracold atoms physics. Moreover, the nontrivial 
dependence of the coefficient from controllable external parameters provides the possibility 
to use these parameters to control the dynamics of the atoms trapped in the optical lattice. 
Thus the key aspect of this paper is the analytical determination of the hopping and interac- 
tion coefficients as a function of experimental parameters such as magnetic trap frequency, 
laser intensity, wavelength, angle between laser sources, s-wave scattering length etc. 

We would like to stress that the procedure followed here for the determination of the 
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coefficients is statistics-independent: the bosonic or fermionic nature of the atoms loaded 
into the trap is completely taken into account by the commutators of raising and lowering 
operators that will be described in the paper. For example with the calculations performed 
here it seems feasible to go beyond the approximations that lead to the Bose-Hubbard model 
in the description of the BEC dynamics in optical lattices, taking into account the specific 
nature of the interaction between alkali atoms in a low density regime. Even for a ground- 
state calculation it can be shown that it is necessary to include levels beyond the single 
particle ground state (see 

The confinement model considered here has a direct experimental relevance (see e.g. 
28|). However, while in a number of atoms of the order of 10^ is considered, 

allowing thus the adoption of a semiclassical model, we focus on a low occupation-number 
regime similarly to what is done in j3| and js^, yet extending to a multi-band model whose 
correctness is limited by the validity limit of the space-mode approximation. 

Challenging tasks for the future will include the determination of tractable yet inter- 
esting models for different aspects for theoretical condensed matter physics and quantum 
mechanics. On the other hand the experimental realisation of systems that exhibit a be- 
haviour which can be described in the framework of the various models here proposed, would 
represent an important achievement for both condensed matter experimentalist and theo- 
reticians: the main difficulties seem to arise form the nearly-single atom trap loading and, 
quite naturally, from the coupling with the external environment. 

Throughout the paper we have tried to emphasise the generality of the procedure followed. 
However we have decided to write down and plot few numerical values of the coefficients 
to stress the fact that this calculation is a direct and relatively simple tool to shape out 
simplified and approximate Hamiltonians for different physical situations. 

In section [H] we depict the potential configuration of the system, moving then to the 
description of our field-theoretical approach. The field operators are written in terms of 
mode raising and lowering operators. Each mode corresponds to a set of quantum numbers, 
one of them identifies the lattice site (hence space-mode approximation) while the others 



describe on-site quantum numbers (local-mode)[3l|. As previously stated this choice is not 
unique, but symmetry constraints suggest expansions that emphasise conservation laws and 
selection rules. 

In section IIIII we evaluate the expression of the Hamiltonian hopping and interaction 
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coefficients and we try to describe the interaction coefficient symmetry properties into some 
detail. 

The purpose of section IIVI is twofold. One the one hand we show how, with suitable 
approximations, the Hamiltonian of the system reduces to known cases, such as the Hubbard 
Hamiltonian or a trivial non-interacting Hamiltonian. On the other we introduce a novel 
Rotational Hubbard Hamiltonian, as a ffist instance of the involvement of higher order 
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approximations. For this case, by means of established mean-field approaches 
suggest a possible path of research involving general group-theoretical procedures [33|. It 
will be shown that these procedures, even if the explicit solution for the ground state is not 
given, allow to grasp interesting aspects of the physics of the model here discussed. 

We have included two Appendices where the relatively simple but lengthy calculations of 
the tunnelling and interaction coefficients are provided explicitly. In Appendix El there are 
various plots of multilevel hopping parameters that supply a good example of the scenario 
that we are moving in and may constitute a good starting point for further investigation. 



II. FERMIONS TRAPPED IN ID OPTICAL LATTICES 
A. General features 

I — I 

The general field-theoretic Hamiltonian (see e.g. 25] ) with 2-body interaction can be 
written as 

H = y"cir^t(r)Hi5(r)^(r) + j drdr'^^{r)i!^{r')H2b{r,r')^{r')^{r) (1) 

where Hif,(r) represents the 1-body term of the Hamiltonian (kinetic -|- external potential 
term) while H26(r,r') the 2-body interaction potential term, \l/(r) is the field operator and 
^^(r) its adjoint. 

As previously mentioned we will stick to neutral fermionic atoms loaded into a ID optical 
lattice. The lattice is generated by two lasers counter-propagating along the x-axis, with 
wavenumber K. The depth - or height, depending if red or blue detuning of the laser is 
considered - in each point x the potential is proportional to the intensity of the laser and 
thus, according to the considered setup, to sin^(2i^'x), for the evaluation of the multiplicative 
constant see e.g. |3|. Here we set the multiplicative constant equal to mu^ / {2K'^) where uj 
represents the harmonic oscillator frequency in the second order expansion of the term V^xt- 

4 



Global confinement is ensured by a cigar-shaped magnetic trap with principal axis along 
the x-direction (see e.g. js^). This trap can be modelled by a 3D harmonic anisotropic trap 
of axial and radial frequencies equal to and Q± respectively {Q^ ^ 

The magneto-optical trap can be thought as if the constituents of the system were trapped 
in the cigar shaped potential with a "slicing" effect of the laser, giving rise to a linear array of 
3D prolate harmonic oscillators. Besides, the radial trapping frequency has a deep influence 
on the interaction among the constituents of the system, allowing to control the volume of 
each "disk". With the previous assumptions Hit becomes 



Hlb — Ekin + Vext 



where 



kin 



2m 

Vext = J + nip'] + ^ sin^Kx) , (2) 

and the second term of Ve^t represents the harmonic confinement of the magnetic trap, while 
the third one corresponds to the optical potential and p' = y'^ + z' . For future convenience, 
we write equation Q as 



Hi, = Ek,n + } Vj+\ v,xt -y y,\ (3) 

with 



o2 -2 " I 2 2 I o2 2 



(4) 



^ji^) — n — j) where n(x) is the rectangle function (n(x) = 1 for — 1 < x < 1, 
n(a;) = elsewhere), k = l±K (with l± = ^Jfij (ma;^)) and Xj = (x— Here the 
harmonic axial confinement of the magnetic field has been considered as a site-dependent - 
with j site index - constant addictive term, merely shifting the local minima of the optical 
potential. From Eq. ^ with the properties of the rectangle function we obtain 

Hi, = J2 nj-(^) [(Ek^n + V,) + {Vext - V,)] (5) 
j 

Hereafter the axial confinement of the magnetic trap will be neglected (small fix)- 
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We are now led to consider two different terms in Eq. The first represents a local 
harmonic-oscillator Hamiltonian 



iho 

j 



2m 2 







(6) 



and a hopping one 



n,(x) 



sm (Ax) X • 



(7) 



27^2 ^ ^ 2 ^ 

The term Vj is the local second-order expansion of the optical potential, thus equation ((Tj) 
represents the discrepancy between an harmonic potential and the true optical potential, 
describing hopping of atoms between neighbouring sites. 

Neutrality of the atoms, ensuring a finite-range interaction allows us to introduce a 
pseudo-potential approximation (see e.g. js^) 

U{r) = y2Uj{x)as6{r)—r, a, := 8 

or m 

j 

where r is the interatomic distance and the s— wave scattering length (a^ in our approx- 
imation is considered constant). The validity of this model is ensured by the low energies 
involved in these interactions, direct consequence of both low temperature limit (virtually 
zero) and diluteness (low Fermi energy). Besides, the form of Eq. (jSJ shows that on-site 
terms only will contribute to the interaction Hamiltonian. Thus Eq. can be rewritten in 
the form 

^ = J2[[ ^^r^^(r)(Hf(r) + Hf""(r)) ^(r) + 

+ ds J drdr'¥{r)¥{r')6{r -r')'^{r')'^{r) . (9) 

B. The (space-|-local)-modes expansion 

The choice of the basis for the expansion of the field operators is crucial. As already 
suggested by the grouping of terms in Eq. ©, we will choose a basis constituted by local 
harmonic oscillator eigenfunctions. In addition, because of the symmetry of the system we 



have chosen central-symmetric 2D h.o. eigenf unctions for the 2D isotropic radial h.o.js^, 
instead of decomposing it in ID h.o. eigenfunctions, this will give us deeper insight into 
conservation laws and selection rules imposed by the symmetries of the system. We then 
have 

^(x) = ^ Un,{x - Xi)Cj^rnip,(p)^i(^)Cn^,J,m,i,a (10) 



i,nx,J,m,cr 



with 



1 - 
Unix) = Hnix/Qe ^ , (11) 

V2"n!V7r/^ 

2im4) / \ 2rri 

CjUp, 0) = ^ [J-J Ly_^ ip/L) , (12) 
Cjm = \J iJ + ra)\ / iJ — m)\ , ^(a) is a spin function and = ^/h/JmUxj- In this de- 



Hn represent the nth Hermite 
with Ly_^ix) a generalised 



composition Unix) is a ID harmonic-oscillator eigenfunction 
polynomial) and Cj^m a 2D harmonic-oscillator eigenfunction 
Laguerre polynomial. 

Fermionic operators will thus have 5 indexes: 3 of them in^, J,m) identify (2+l)D local 
harmonic oscillator states, while i identifies the site and cr the spin. While has its 
usual interpretation of ID harmonic oscillator number operator eigenvalue, J and m can 
be construed as angular momentum and x— axis component of the angular momentum, 
respectively. 

This decomposition can be thought as a generalised space-mode approximation, with 
additional local modes that, in the present case, correspond to the local (2+l)D harmonic- 
oscillators quantum numbers. If not explicitly required, we will use (pa = Un^ix — 
Xi^)Cj^^niaiPi'P) ii^a)i with a = {ua, Ja,™a,'ia,o'a} to simplify the index notation. We 
wish to stress that decomposition is an approximation of field ^(x): there is a non-nil 
overlapping between wavefunctions belonging to different sites, thus orthogonality is not 
fulfilled. Nevertheless these overlapping integrals are supposed to be small, ensuring the 
consistency of this choice [lol |. 

In the forthcoming calculation of the interaction term, Eq. (jlip allows us to easily 
recognise that m is a conserved quantity. If we come back to ((H), with the decomposition 
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(fTUI) we obtain 



dr<p:{r)H';"Mr)cicp+ / t/r 0:(r)Hf->^(r)ctc^+ 



7,(5 



(13) 



III. HAMILTONIAN COEFFICIENTS 



We are now in the position to calculate all the coefficients in Hamiltonian (|13|) . The ffist 
term becomes 

/oo 
cirn,(x)C(r)0/3(r)ctc^, (14) 
„ „, „ -oo 



where A/j is the (2+l)D harmonic-oscillator eigenvalue 



A 



nfs,Jfs,mi3,ii3,af3 



nio,[np + - ] + h{2Jf^+ 1) 



Eq. p4j] can be written as 



j,a,f3 a 



rir 



(15) 



(16) 



where the second Kronecker delta is a consequence of the space-mode approximation, i.e. we 
consider only superposition of wavefunctions among which at least one is a local harmonic- 
oscillator eigenf unction, while the ffist one stems from the orthogonality of the (f)-y{x) func- 
tions. 

We move now to the evaluation of the integral in the second term of equation (|13|). 
Namely 

f^tunn = ^ I dr0:(r)n,(x)Hf""(x)0^(r) (17) 

being H*"""(x) independent of radial and spin degrees of freedom, we can rewrite equation 
(|T7i) as 



H 



tunn 



[x)Unu,iij {x)ci^Cf3. 



(18) 



With the same assumptions of the local harmonic-oscillator case the integral in equation 
f|TK|l becomes 



Kn^,n,nw^ / dye-^H,,^ (y - r) Hf '^"(y)e--if„, (y) 



(19) 
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where Kn^ 



,n0 



[2'^"~^"-i^nJ it] -^/^ and we have put y = {x — i^d)/!^ (where d = n/k), 



T = {ip — ia) and = Kl^. By substituting the expression of H*"""(?/) from equation Eq. 
((Tj) we obtain 

Kn^^np^x j dye' 
If we define 



Hn^ {y - t) Hn, {y) 



1 - cos(2fiy) y2 



c%- (20) 



Ta,l3 — ^;^Sj^Jij6m^^rnf,Sa^,a0fii^x / dye' 



X Hn, (y) Hn^ {y - t) 



y^ 1 - cos(2f^?/) 



(21) 



equation (fT^ becomes 



H*""" = -)^T.,/.at,£/3. (22) 
Then the term Tq, c^cJ^Cq can be incorporated into the \\^-°- term, giving 



/Iq, Aq, XJ;^ 



(23) 



We will here skip the explicit solution of the integral in Eq. ()21|). along with the analytic 
expression of T, which can be found in Appendix^ These calculations allow us to write 



-^a,/3 ^JctjJf) ,mp ,a-p -^ric ,np ,ia ,ii3 ■ 



(24) 



In fig. Qwe have the plot the coefficient Tn^^ni3,ia,ii3 as a function of the ratio between distance 
and the period of the optical lattice, for nc, n^j = 0, 1. In boldface we have marked the points 
corresponding to discrete values of the ratio x/d, i.e. the points with a relevant physical 
meaning. The values of T plotted here are in arbitrary units. Even if the correctness of 
the above procedure seems undoubted, it must be remembered that it is entirely based on 
the space-mode approximation, whose validity depends on the overlapping of wavefunctions 
belonging to different sites and thus might be violated. 

These plots show how the tunnelling amplitude varies with the distance. In particular it 
is clear how, for long-distance tunnelling, there is a negative exponential dependence. Never- 
theless, if the experimental conditions are properly chosen (i.e. angle between counterprop- 
agating laser beams and their power), it is possible to obtain conditions where, for instance 
nearest-neighbour and next-to-nearest neighbour tunnelling coefficients have opposite signs 
(see e.g. Fig|2lTo,o,i^,i^), and thus the model, in that case, might exhibit frustration. 






FIG. 1: Plots of the tunnelling coefficients from Tqo to T22. A detailed discussion of the analytic 
expression of the hopping parameter will be given in Appendix 1X1 

We will now move to the determination of the interaction term, namely the last term of 
equation (fT^ . As a first step, we can write the integral in cylindrical coordinates 



dsj rfrrfr'C(r)0;(r')5(r-r')0/3(r)05(r') = 

j" dxdx u^^i^x x)it^^(x Xj^)it,^^(x -^j^) 

5{(t)- (t)')Cji,^mp{p, (t))Cj^^rn,{p\(t)') 



X / dpdp — 
vr 



(25) 



with p = p/lp and the identity 



5{v) 



5{p)m 



Tip 



As we are dealing with a short range interaction modelled by a 5(r — r') function, we will 
consider on-site interaction only (xj^ = = Xi^ = Xi^). 

This choice is completely justified because the interaction term is modelled by a pseu- 
dopotential term for which nearest-neighbours interactions become negligible. In this case 
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the first integral on the left-hand side of Eq. ()25|) becomes 



1 /2-("Q+»^/3+"7+"i) f , 

U^ = — J ■ — , , , / dxHn ix)HnAx)Hn (x)HnAx)e-^'' (26) 

with X = x/lx, whose explicit calculation is given in Appendix^ Here we just give the final 
result 

(27) 



i^ii - m) I ^ 



with 



,0 if odd 

Eis) = { 28 
UeM7s)HsM Use even 

The summation is to be intended as 4 separate summations over the components of a vector 
s = {sa, S/3, s^, Ss} from {0,0,0,0} to n = {n^, n^, n^, n^}. The norm ||x|| is a 1-norm 
(ll^ll = J2e G = (y-^ (3,1,5) and the 5 in Eq. (PBjl represents the parity selection rule, 
obtained from the explicit calculation of the integral. 
For the radial part of the integral we have 

^P = j j '^P'^^ ^^X,m^ (P' ^) ^\.ra^ (P' (P, 0) I^J,.rn, (P, 0) (29) 

with the definition given by Eq. (jllj) . we can easily perform the angular integration and we 
obtain 

2(5 /""^ 
= / dppRx,^^{p) RX^^^ip) Rj^mM Rjs,mAP)- (30) 

^ Jo 

The reader is again addressed to Appendix ^ for the explicit evaluation of the integral in 
Eq. (jHUI) . The result is given by 

^ ^('^'"''^) 211^11+3/2 (31) 

5=|m| 

with 

fc^j;,,* (^<>-*)!(» + n>ii)! (*- ">»)! 
and, following previous notation, we obtain q = {qa,(li3,<l'r,(ls}, J = {Ja, J/s, J-y, Js} and 
|m| = {|mo,|, |m^|, Im^l}. The overall interaction coefficient can then be written as the 
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product of Eqs. dHD) and (gZj) 



jj . as ^ {J,rn,q) S (g) 



q=\m\ 



V2' 



xr ||g 



l|5||+2||g|| 



+ 11 



(33) 



We are thus enabled to rewrite Hamiltonian (|T^ in terms of the calculated coefficients 
obtaining 



(34) 



_ Ci a, 13 a,/3,7,(5 

We will refer to ()34|) as the generalised Hubbard Hamiltonian. For sake of simplicity, in 
Eq. we have not written down explicitly the selection rules imposed by symmetry 

constraints (see below). 



A. Symmetry properties of the interaction term 

In addition to global symmetry properties, such as 1) rotational symmetry along the 
X-axis and 2) left-right symmetry, reflected by momentum x-component conservation and 
parity conservation for the ID harmonic oscillators along the x-axis, it is clear from equation 
flHHj) that the coefficient Ua,(3,'y,& has some symmetry properties: a) U does not depend on 
the sign of with x = 7,^, provided the conservation m during the interaction 
{rria + = + m^); b) U possesses a permutational symmetry, namely 

We would like to draw reader's attention to the two 6 functions in equation (j33p which 
make explicit the conservation laws that might have been expected by simply considering 
the symmetry of the problem. The first one represents parity conservation, while the second 
conservation of the x component of the angular momentum. In table HI we give the analytical 
value of U for interaction between particles belonging to the first three shells of the 2D radial 
harmonic oscillator and to the first level for the axial harmonic oscillator. These symmetry 
constraints allow to class the possible quantum numbers of the interacting particles according 
to the value of the corresponding U. For example for 

a = {0,1,0, i, a}, /3 = {0, 0, 0, a'} , 
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(5 = {0,0,0,z,cr}, 7 = {0,0, 0,2, a'}, 

and 

a = {0,0,0,i,a} , /3 = {0, 0, 0, i, a'} , 

S = {0,l,0,t,a}, 7 = {0,0, 0,2, a'}, 

we have the same value of U, henceforth the class definition of table We would like 
to point out two aspects of this example. First of all it may be noticed that the angular 
momentum J is not conserved throughout the interaction: this is a general feature of the 
system considered, there is not global rotational symmetry but only in the plane orthogonal 
to the ID optical lattice. Moreover in this particular interaction the value for the coefficient 
U is negative, this appears to be a rare (but not unique) situation. The implications of this 
condition will be pointed out in section IIVI 



IV. SPECIAL CASES 



In this section we derive three model Hamiltonians for fermions in optical lattices. We 
consider, for both cases, only the lowest-state axial quantum number (i.e. = = 0). 
Hence T^^^ can be written as 

Ta,l3 = ^Jc„Jf}^ma,mf,Saa,ai3Tofi,i^,ii3 (36) 

As far as an ultracold gas is considered, it seems feasible to restrict our analysis to the first 
few levels above the ground state (i.e. Jq, = 0, 1/2, . . .). As a first example, we consider the 
case having = as the only radial level allowed and the fermionic gas is spin unpolarised. 
From Eq. (jSSI), along with the previous assumptions, we obtain [^31 

with T = ro,o,jc<,ja+i which is easily recognised as the Hubbard Hamiltonian, whose role in 



the ultracold atoms physics has been pointed out elsewhere Note that in this example 

we have made the assumption that the tunnelling coefficient is significantly different from 
zero only for nearest-neighbouring sites. Nevertheless more involved situations may arise, 
suggesting interesting physical features, as it is shown in Appendix 1X1 
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TABLE I: Values of Ua,i3,'y,5 = [h'^CLs/ {mlxTT'^l'j_)] ^ Ua,i3,'y,s for n/3, n^, n^} = {0,0,0,0}, ia = 
il3, a and a' satisfy symmetry constraints. The value represents the equivalence class described 
in the text. 
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If we now consider a spin-polarised gas in a (radial) multi-level system we obtain 

H=Y1 t^n,inn,i -TY^ {ct/nA+l + 4,i+lCn,i) (38) 



where the absence of the interaction term is related to the symmetry properties of the 
coefficient Uap-yS- The Hamiltonian (pH|l is readily diagonalised to yield 



n n,k,a 



(39) 



14 



with the same procedure followed in the strong coupling limit in the Hubbard Hamiltonian. 
Rotational Hubbard Hamiltonian 

As the last example, we derive a third Hamiltonian that may give the reader a first 
insight on the increasing complexity if higher single-particle levels are taken into account. 
Here we consider a situation where we allow = 0,1/2 (but always = 0, leading to 

Tn^,ia,,n^,i^+i = T as already pointed out) 

1 

H2-level = ^ ^ I^Aarij^a^o- — T ^c|_^ ,^Ci+i,a,<T + 4+l,a,cTCj,a,o- j j + 
i,a a=—l 

Ua,b,c4Y^i^^<^^ib,a'Ci,c,a'Ci,d,a (40) 
i a,b,c,d a, a' 

where Ua^b,c,d = Ua,i^;h,ifi;c,i~,;d,is- The label a (as well as 6, c, and d) has been introduced 
to represent the triplet of harmonic-oscillator numbers {ua, JajiTT-a), hence a = {a,ia,cra}- 
One should recall that originally a = [ua, Ja,fna,ia,o'a)- Here, however, it is convenient to 
write in an explicit way both spin indices a^'s and site indices i^s. 

The triplet a = {ua, JaynT-a) is such that the value a = corresponds to (0,0,0), a = 
1 (0, 1/2, +1/2) and a = —1 —>■ (0, 1/2, —1/2). The axial quantum number Ux has been 
"frozen" to due to the disk-shaped potential form (i.e. 3> u!±) while the radial quantum 
number J has been limited to the values {0, 1/2} as a first approximation beyond the J = 
(Hubbard Hamiltonian see ()37|)). The present model thus enriches the dynamical scenario 
by introducing modes that takes into account the simplest possible rotational processes for 
fermions confined in a well. 

The wealth of the scenario depicted in Eq. ()4()|1 arises from the level-dependence of the 
interaction coefficient Ua^b,c,d- In fact Uafi,c,di as a function of the energy levels may provide a 
useful tool to simplify Eq. ()40|) hinting the best strategy for both numerical and analytical 
analysis of this model. Two main aspects concerning these coefficients are worth repeating 
here: a) the m^-conserving nature of the interaction, related to the symmetry properties 
of the confining potential and of the interaction coefficient (see section UlI A|l . reduces the 
number of possible processes; b) the symmetry properties of Ua,p,-^,& (see Eq. fjHH|) ) allow the 
grouping of interaction terms, accordingly to what has been done in table HI 

From a general point of view, in Hamiltonian ()40|) the hopping factor may be construed 
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as a multichannel tunnelling coefficient, where the radial quantum numbers identify the 
channel label, in the same spirit of Hamiltonian (jHIIj) . Incidentally, this is true if axial 
degrees of freedom are "frozen" to = = 0, otherwise there is tunnelling among levels 
with ria 7^ Up, for some a and (3. 

Hamiltonian (j4Up can represent a situation where single traps are loaded with a small 
number of atoms, as to fill the first two radial levels of the local harmonic oscillator. To 
experimentally obtain one of the different simplified Hamiltonians - like ()40|) - it is necessary 
to have control of four experimental parameters: laser intensity, angle between counterprop- 
agating lasers, axial magnetic trapping frequency, scattering length. With these parameters 
it is possible to gain full knowledge of "lattice constant", interaction parameter, shape and 
depth of the 3D harmonic traps. The most critical point seems the few-atoms loading of 
the trap but a technique involving a 3D anisotropic array -a sort of 2D array of ID arrays- 
might overcome the problem. 

In this picture, the interaction coefficient U can then be used as a source of entanglement 
between different channels. Moreover the possibility of experimental control of the scattering 
length, and thus of the interaction term, via an applied magnetic field may provide an useful 
tool of external manipulation of the state of the system in the rich scenario here depicted. 

To outline future paths of research, we will here sketch a way to set up a mean-field 
procedure for the Rotational Hubbard Hamiltonian. The main interest of this approach 
resides in the possibility of a general discussion of some features of the model which have 
a direct experimental relevance. For example it is possible to state that, according to what 
is usually affirmed in the literature , no BCS-like ground-state is possible for repulsive 
interaction, while for an attractive two-body potential a paired ground state is possible. The 
flexibility of experimental techniques involved in the study of ultracold atom physics allows 
to envisage experimental conditions where these two different regimes are attained. For 
example exploiting a Feshbach resonance it is possible to drive the scattering length from 
positive to negative values leading thus the system through a quantum phase transition. 

The analytic procedure adopted hereafter deeply relies on the concept of quasi-free state 
3^ . In our situation the following definition of quasi-free state can be adopted: 

I) all correlation functions can be computed from Wick's theorem; 
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II) four fermionic expectation values over a quasi-free state have the form 

< 016162636410 >=< 0|6ie2|0 >< 01636410 > - 

< 0|6ie3|0 >< 0I6264I0 > + < 0|6ie4|0 >< 0I6263I0 > 



with Cj = Cj, c|. In particular we would like point out how the three terms on the right 
hand-side will lead to the direct, the exchange and pairing energy term of a Hatree-Fock- 
Bogoliubov mean field Hamiltonian, which, for our RHH becomes 



fjHFB 
^2-level 



a,b,cA 



i,a,b,c,d, 



with 



i,a,a 

Xi,a,a,b,a' =< 4>HFB\cla^^Ci^by\(f)HF > 
^i,a,a,b,a' =< (pHFB\Ci^a,aCi,b,a'\4>HFB > 

The set of generators — ^Sajsi^ < 0: ^ f3 < r), CaCjj, c^c'i^il < a 7^ < r) 

the following commutation relations 



obeys 









- 2^H^ 




5jk{.c\ci - 


-5,7 

2 ti 










- ^AA 




Ci Cj , 






5ik{ 






) + 





kj{c{ci - l/25ki) - 5ii{c{cj - l/25kj) 
5ki{c\ci - 1/24) 



(41) 



allowing to state that the dynamical algebra of this new Hamiltonian, which is now quadratic 
in terms of q, c| can be easily recognised to be so(2r) js^. 
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Having determined the dynamical algebra of the model Hamiltonian, enables us - at least 
in principle- to find the ground state of the system with a straightforward procedure. As 
it will be clear form the subsequent discussion, the main difficulties arise as the number 
the generators of the so{2r) algebra grows with r(2r — 1). For instance for the two-site, 
J = 0, 1/2 model, the Hamiltonian dynamical algebra will have 276 generators. 

In spite of the technical difficulties (both analytical and numerical), it is appropriate to 
apply algebraic techniques to diagonalise H2-^^ei- -^^ stated before, this general approach 
will give some insight to the ground state properties of the system. If we consider a unitary 
transformation g G S0{2r) we can write 

Hd = 9HHFB9-' (42) 

where Ha is diagonal. As a direct consequence the ground state \(j)HF > of H2Ji^ei t)e 
written as 

\(f>HF>= g\0, 0,0,..., 0,0 >=g\0> (43) 

where |0 > can be defined as the Bogoliubov particle vacuum (ground state of Hd). Following 
> represents a possible choice for the extremal state for the S0{2r) group with U{r) 



3, 



as the corresponding maximum stability subgroup. Leading to 

^|0 >= nh\0 > n\0 > e'^^'^'^ (44) 

where: 

exp ^ (ria,pclcl - H.c^j G (45) 



n 

The phase appearing in Eq.()44j) has no relevance for our purposes, as we are interested in 
the evaluation of observable expectation values. 

The problem mentioned above about the size of the dynamical algebra, appears here with 
all its implications. It is necessary to exponentiate the operator X]i<a^/3<r {va,pc^a^f^ ~ H.c^ 
which is a 2r X 2r matrix in the faithful matrix representation. 

Nevertheless, for a repulsive two-body potential, the pairing term can be neglected j^^l, 
thus the dynamical algebra of the system becomes U{r). Following Q], we can express the 
Hamiltonian ground state as 



')HF >= exp ^ {r]a,f3cicf3 - H.C.) \0 > (46) 

A:+l<a<r 
l<j<k 
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where 

|0>= |M^_^,0,...,0> (47) 

k 

which is, in fact, the ground state of the non interacting Hamiltonian. It is worth noticing 
that this general procedure can be greatly simplified if further constraints, related to sym- 
metries of the problem, are imposed onto the coefficients rja^is- For example, if we consider 
the two-site(y4,i?), J = 0,1/2 case, due to equation (jHBj) the matrix f] with elements rja^^ 
will have the form 






^71,2 


^1,3 


Via 


^1,5 


Vi,e 


^71,7 






















^2,3 


V2A 


V2,5 


V2,6 





^72,8 














-111,3 


-^72,3 





V3A 


113,5 


1l3,6 








^3,9 











-Via 


-^2,4 


-V3A 





Vi,5 


11^,6 











^4,10 










-^2,5 




-^4,5 





115,6 














^75,11 







-V2,6 


-V3,6 


-^4,6 


-V5,6 




















^76,12 








































-V2,8 






































-^3,9 






































-^4,10 






































-'75,11 























V 














-^6,12 





















(48) 



thus impressively reducing the computational effort needed to evaluate the exponential in 
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equation ()45|) . In Eq. ()48|) we have have assumed the following convention 

{n„ = 0, Ja = 0, rria = 0,ia = A, aa =T} ^ 1 

{n^ = 0, Ja = 0, TUa = 0,ia = A, =[} ^ 2 

= 0, J„ = 0, m„ = 0, z„ = 5, a„ =t} 7 
{n„ = 0, J„ = 0, m„ = 0, 2„ = 5, da =1} 8 



V. CONCLUSIONS 

In this paper we have investigated the complex structure of fermion interactions for a 
fermion gas distributed in a linear periodic array of potential wells. Based on the stan- 
dard many-fermion quantum field theory endowed with a potential distribution mimicking 
a realistic experimental setup, we have calculated analytically the hopping and interaction 
coefficients that describe the interactions of fermions within a generalised multimode Hub- 
bard Hamiltonian. Their dependence on the external controllable parameters (such as laser 
intensity, magnetic trap frequency, wavelength, and scattering length) has been determined. 

Our analysis shows that, except for two particularly simple cases (the gas of spin unpo- 
larised fermions and the gas of noninteracting spin polarised fermions) , models with different 
degree of complexity can be derived depending on the interaction processes one decides to 
account for or to neglect [consider, e. g., that, in principle, one might introduce an unlim- 
ited number of (local) rotational levels]. In this respect, our simplest nontrivial model ()4()j) . 
which is able to account for the (local) rotational activity of fermions, appears to be far 
more complex than the Hubbard model or the spin-polarised noninteracting model derived 
in section HVl 

Therefore, the first objective of our future work is to perform a systematic study of 
model ()4()j) . Based on the present analysis and exploiting the interaction-parameter scenario 
here depicted, the second objective is to recognise the significant regimes characterising the 
confined fermion gas and to derive the relevant models from Eq. fl34|) . 

We would like to stress once again how the analytical knowledge of the coefficients in 
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principle allows us to tailor Hamiltonians performing specific tasks. 

An aspect that certainly deserves attention is the study of the zero-temperature phase 
diagram of model (j4Up (and, more in general, of sufficiently simple -and thus tractable- 
models derived from the gHH) and of the relevant phenomenology aimed at suggesting new 
possible experiments. To achieve a reliable description of these systems, several established 



analytical and numerical approaches (see e.g. j38|, |39|, and 



mm 



42j . respectively) can be 



4lL l42il . Moreover in the recent 



implemented in analogy to what has been done for bosons 

past several authors (see e.g. [4^ ^|) have proposed to use entanglement measures as a 
quantum phase transition identifier. We think that our model can represent a good test-field 
for this new approach to quantum-phase transitions. 



APPENDIX A: TUNNELLING COEFFICIENT CALCULATION 



In the following calculation we will fix np > Ua, without loss of generality, as it can be 
easily verified. 

The integral in Eq. !^21\i can be decomposed in the sum of three terms 

Qn.,., ^ Q..,n, ^ ^ J ^y^-^^^^ _ ^) /(^)e-^i7„^ (y) (Al) 

with 

^ f l - cos(2r]y) _ y^' 
[ 41]2 2J ' 

integral ()A1|) becomes 



er'"" = / ^^""^^ co«(2fi2/) X (y) Hn, (y-r), (A3) 

er'"" = -ll dyy' e-^^^ X H^^ (y) H^^ {y - r) , (A4) 
The substitution ( = y — t/2 yields 

er =C^nJ dCe-^'H^^ (c + ^) H^, (c - ^) (A5) 
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where = e '^^/^/(4f2^). We then use the Hermite polynomial identity 

H^{x + y) = J2(f)Hk{x){2yr 

to obtain 



(A6) 



e 4 
4lF 



Z,fc=0 



With the orthogonality of Hermite polynomials 

) 

we are able to perform the C integration 



(A8) 



-2)'/!. 



(A9) 



It is worth noting that the summation extends to n^-, that is (^) = if a < h. From 
it can be verified that the last summation is related to generalised Laguerre polynomials, 
giving 



We move now to the calculation of O2"' which is given by 



e 



1 



dyexp 



r 
2 



(y - rf 



H^Ay-r) cosily). (All) 



Hn^ (y) exp 

Eq. ()A11|) . with the substitution ( = y — t/2, can be written as 

er = exp (-rV4) J dCe-<'H^^ (c + ^ H^, (c + cos m . (A12) 

Again, using Eq. ()A6|) gives 



e 



E 



Lk=0 



x(-l) 



na — k 



dCe-^"Hu{0Hi{C)cos{2^C) (A13) 



The integral in equation ()A13|) can be interpreted as the real Fourier transform of the 
function 
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Recalling that 

:F[f{x)g{x)]=:F[f{x)]*r[g{x)] , 
where JF [■] indicates Fourier transform and * convolution product, we obtain 



e 4 I rir 



41]2 ^ V / 

Lk=Q 



k 



X 



-1) 



T[e-'-Hk{0]*J'[e~'-m{()\ 



(AM) 



giving 



4^2 ^ \ I 

Lk=0 



X Me 



dee 



. (A15) 



-i^fc(e)i^Ke-2^^) 

The integral on the left-hand side of equation ()A14|) can be solved following the same pro- 
cedure used for ©""'"''^ 

j dee-'''^Hk{e)e~'^'-^''^''^Hi (e - 2fi) = (-l)^'-'v^e-^ {2nf-' /!2'Lf-' (41]2) (A16) 



giving 



e 4 



E (7 "-(-!)- 



X 



The calculation of ©^"'"''^ is quite straightforward. Applying twice the identity 
we can write 0^"'"'' as 



(A17) 



(A18) 



e 4 



dC^e 



1 / r\ . 2n„ + l 



X^n.-2(C+0]^n,(C-0 . (A19) 



With the same procedure used for 9""'"'''^ Q^t*."-/? given by 



^2, 



v/^2>„!r"''-""e-"r 



K + l)K + 2) 



2 -^nc«+2 



ni3-na-2, 2 



(r72) + 



+ 4 ^n.-2 I 2 i + 2 



1 j-np-na I T 



. (A20) 
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Hence T^^p becomes 




(A21) 



— 



with e""'"^ e2"'''^ and 63"'"'' given by formulas (IMo|) . (IXTTI) . and (IMol) . respectively. 

We have here the plot of T„^,n^,i„,ic«+i &s a function of the difference ia — ip for values of 
Ua and ranging from to 2. 

The long distance exponential decay is common to all tunnelling coefficients, regardless 
of the energy level. On the other hand its detailed shape has deep relevance for nearest- 
neighbours and next-to-nearest-neighbours (i.e. there may be sign changes passing from 
Tn,m,i,i+i and T„ m,i,i+2) as shown in figures (0111). Another interesting feature is that an 
extra-term, due to "on site" tunnelling coefficients, bust be added to the harmonic-oscillator 
energy term. 



FIG. 2: Plot of Tn^^np from To.o,j^,j^ to To^2,ia,ifj- The solid line represents the case of a ground-state 
tunnelling. In this case the hopping parameter To^o,ia,i/3 is always positive. However, inter-level 
tunnelling already shows sign changes. 

APPENDIX B: INTERACTION COEFFICIENTS 

We provide here the detailed calculation for the interaction term matrix elements. To 
solve integral (PU)) 




X 
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FIG. 3: Plot of Tn^^np from Ti^o,ja,j^ t° 2^i,2,jQ,j^- this situation the intra-level tunnehing term 
^i,i,«Q,j^ is always negative, but if a different external parameter choice is considered, the sign 
change can be placed between 1 and 2. 




FIG. 4: Plot of Tr 



a, rip 



from 



2,0,ia,ip 



to 



. In this case the sign change for the intra-level 



tunnelling term occurs for the specific parameter choice performed here, but it can be removed by 
a different choice of the external parameters. 



we exploit again Eq. ()A6|) . obtaining 



' 2-("-a+"/3+"7+"«) 



E 



71/3 \ I 



na+rifs+nj+ns-iia+ifi+i-i+is) ^-'^('^ 



fBll 
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In previous equation the summation must be intended over four independent of sg = . . .ng 
with 9 = a, P, 7, 6. With the substitution 



I daKTe-'^" = 5«,2NV^°^'r [{a + l)/2] 



(B2) 



('^a,2N indicates that a must be an even number) Eq. ()B1|) becomes 



E 



V2' 



]|5||+3 



r 



\n\ 



lln||,2N 



with 



and 



= {oq., a/3, a^, a^} 1-norm: ||a|| = 



ag 



n 



1 /ne 



HsM 



(B3) 



(B4) 



(B5) 



where S (s) = for odd ig. The 5 function in Eq. ()B3|) should be written as 5{||n||-||5||,2N- 
However the condition \\ig\\ =even aheady imphes =even. We are then allowed to write 
in Eq. (IE3|) : 

^||n||-||s||,2N = <^||n||,2N • 

We solve now the radial part of the interaction term integral written in Eq. (j29|) 



where rj = {p,4>) and d'^r] = pdpdcj). Following [36[, we express Cj^^mSPi 'P) i'^ terms of a 
finite sum: 



1 



X 



2go 



qa = \m,a\ 



(B6) 



Hence the radial part of Cj^^maiP: 4>) can be written as 



ga = |mai| 



where 



A. 



(B7) 



(B8) 
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Substituting Eq. (jB7|) into Eq. (j29I) we have 



25. 



Ja J/l J~1 Js 



qa = \ma\qi3 = \mi3\q-, = \m-^\qs=\ms\ 



dp:^[j-J (B9) 

that, with the same notation of Eq. ()B3jl . becomes 



^^2 2^A(J,m,gj ^llgll+a/^ ^ 



g=|m| 

with 



A(J,m,g)= J] A,. (Bll) 
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